arXiv:1503.02694vl [astro-ph.EP] 9 Mar 2015 


Draft version March 11, 2015 

Preprint typeset using D-T^X style emulateapj v. 5/2/11 


ARE PROTOPLANETARY DISKS BORN WITH VORTICES? - ROSSBY WAVE INSTABILITY DRIVEN BY 

PROTOSTELLAR INFALL 

Jaehan Bae 1 , Lee Hartmann 1 , Zhaohuan Zhu 2,3 
Draft version March 11, 2015 

ABSTRACT 

We carry out two-fluid, two-dimensional global hydrodynamic simulations to test whether proto- 
stellar infall can trigger Rossby wave instability (RWI) in protoplanetry disks. Our results show that 
infall can trigger the RWI and generate vortices near the outer edge of the mass landing on the disk 
(i.e. centrifugal radius). We find that the RWI is triggered under a variety of conditions, although the 
details depend on the disk parameters and the infall pattern. The common key feature of triggering 
the RWI is the steep radial gradient of the azimuthal velocity induced by the local increase in density 
at the outer edge of the infall region. Vortices form when the instability enters the nonlinear regime. 

In our standard model where self-gravity is neglected, vortices merge together to a single vortex within 
~ 20 local orbital times, and the merged vortex survives for the remaining duration of the calculation 
(> 170 local orbital times). The vortex takes part in outward angular momentum transport, with a 
Reynolds stress of < 10 -2 . Our two-fluid calculations show that vortices efficiently trap dust particles 
with stopping times of the order of the orbital time, locally enhancing the dust to gas ratio for particles 
of the appropriate size by a factor of ~ 40 in our standard model. When self-gravity is considered, 
however, vortices tend to be impeded from merging and may eventually dissipate. We conclude it 
may well have that protoplanetary disks have favorable conditions for vortex formation during the 
protostellar infall phase, which might enhance early planetary core formation. 

Subject headings: accretion disks, hydrodynamics, instabilities, stars: formation, stars: pre-main se¬ 
quence, waves 


1 . INTRODUCTION 

Dust particles in protoplanetary disks feel drag forces 
from the gas. In general, the disk gas has out¬ 
ward pressure support and so rotates at sub-Keplerian 
speeds, whereas grains attempt to move in Keplerian 
motion. This effect is most significant for the parti¬ 
cles that have the stopping time comparable to the or¬ 
bital time, t s ~ 1 , where t s and fI denotes the stop¬ 
ping time and the dis k rotation frequency, respectively 
(jWeidenschillind fl977fl . Radial drift of dust is thus de¬ 
pendent on the disk conditions under which the parti¬ 
cles reside in, but for centimeter to meter-sized parti¬ 
cles the resulting radial drif t velocity usually reaches a 
few tens to hundreds m s -1 (|Weidenschillinelll977l Il993l : 
iKlahr fc Bodenheimeri 120061 ). The corresponding radial 
drift time at 10 AU, for instance, is less than a few 
x 10 3 years which is much shorter than planet formation 
timescale as well as disk lifetime. 


Formation of vortices in protoplanetary disks could 
be important in the context of planet forma¬ 
tion due to their ability to efficiently trap dust 


particles (Adams & Watkins 1995t Barge & Sommeria 

1995 

: iTansia et alJ 1996j; Bracco et al.l 1999 r |Chavanis 

20001 

: Fromang & Nelson 20051: Inaba & Barge 2006: 

Hcng & Kenyon! 201(1 

Birnstiel et al.1 20131: IZhu & Stones 

20ll Zhu et ah 2014 

). Anticyclonic vortices - with 


highest pressure at their centers - are especially irnpor- 


jaehbae@umich.edu, lhartm@umich.edu, zhuzh@astro.princeton.edu 

1 Deptartment of Astronomy, University of Michigan, 1085 S. 
University Ave., Ann Arbor, MI 48109, USA 

2 Department of Astrophysical Sciences, Princeton University, 

4 Ivy Lane, Peyton Hall, Princeton, NJ 08544, USA 

3 Hubble Fellow. 


tant since they can survive long (hund reds o f orbits), 
while cyclonic vortices dissipate quickly dGodon fe Liviol 
ri999h . Also, at the core of anticyclones, the radial and 
azimuthal gas pressure gradient s vanis h and the refore 
gas t here is in Keplerian motion (IKlahr fc Bodenheimeri 
12009 ). This is an interesting feature of anticyclones be¬ 
cause they not only concentrate dust, but once parti¬ 
cles reach to the vortex center there is no drag force if 
turbulence and scattering of particles is absent. There¬ 
fore, vortices can be regarded as promising sites of 
planet formation. In fact, highly non-axisymmetric, 
vortex-like structures in protoplanetary disks have been 
observed in submillimeter and millimeter interferomet¬ 


ric observations: LkHa 330, SR 21N, HD 135344B 


(Brown et all 2009 

), Oph IRS 48 (van der Marel et al. 

2013 

), HD 142527 

Casassus et al. 2013c Fukagawa et al. 

2013 

), SAO 206462, SR 21 (Perez et all 20117 and 


PDS 70 (iHashimoto et al.ll2015f l. 


Possible ways to form vortices in protoplane¬ 
tary d isks are throu g h th e Rossby wave i nsta bility 
(RWI; ILovelace et ahl HOOl ILi et, all l2000l l200l or 
the Papaloi z ou-Pring l e inst ability (|Papaloizou fe Pringlel 
11984 119851 lHawlevI I1987T) . Generally, the instabili¬ 
ties are triggered by the additional shear arisen from 
the change of azimuthal velocity gradient over the ra¬ 
dius of interest. In two dimensions, it has been both 
analytically and numerically shown that the RWI can 
be triggered at a minimum of the generalized vorten- 
sity ti = k 2 /(2YUS' 2 / 7 ) (|Lovelace et all 119991 ILi et al.1 
2000 ), where k is the epicyclic frequency, T, is the sur¬ 
face density, Q is the rotation frequency, and S is the 
entropy with the adiabatic index 7 . Recent numeri¬ 
cal simulations show that the RWI appears to trigger 
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in a simil ar way in three dimensions, as in two di¬ 
mens i ons (iMeheut et al.1 [20101 l 2012 at ILvra fc Mac Lowl 
l2012t iRichard et al.l 120131 : iLinl 120141) . The question is 
then how to make a vortensity minimum in protoplan¬ 
etary disks? Previously, it has been suggested that 
a vortensity minimum can form and the RWI is trig¬ 
gered at the edge of the disk d e ad-zone (llnaba fc Bared 
20061: lYarniere fc Tagged 120061: ILvra et all I2008L 120091 : 
Lyra fe Mac Lowll2012t lRegalv et al .1120121) or at the edge 

o f ga ps carved by a planet (jde Val-Borro et al.ll2007i : iLinl 

1ml . 

In this paper, we propose another mechanism that 
possibly drives the RWI in protoplanetary disks: pro- 
tostellar infall. This is partly motivated by the recent 
ALMA image of HL Taro showing multiple ring and gap 
structures which might be generated by planets. We fo¬ 
cus on the fact that HL Tau is embed ded in an enve¬ 
lope f rom which it still acc retes material (IBeckwith et al.l 
119891 : lHavashi et al.lfl993l) . If the gaps in the ALMA im¬ 
age are indeed a signature of planets, this implies that 
planet formation can happen in the very early phase of 
protostellar evolution. Our simulations show that proto- 
stellar infall from natal cloud can generate local vorten¬ 
sity minimum near the centrifugal radius, inside of which 
radius infalling material falls onto the disk, and can trig¬ 
ger the RWI. While the RWI activity depends on the 
characteristics of infall model (e.g. radial infall profile, 
existence of shear between infall and disk material) as 
well as disk parameters (e.g. viscosity parameter), we 
find that the RWI can be triggered under a broad cir¬ 
cumstance and is a possible way to form vortices. Our 
results suggest that vortices can form during very early 
evolution of protostellar systems, and this may ease the 
timescale problem for giant planet formation. 

2 . NUMERICAL METHODS 
2.1. Infall Model 

The infall model is implemented by adding the corre¬ 
sponding terms to the hydrodynamic equations as below. 

r)V 

+ V • (S g u s ) = E in (1) 

^ 9 + v 9 ' ^ v g^j = — SgV^+V-Hg+Ein (2) 

In the above equations T, g is the gas surface density, v g 
is the gas velocity, P g = E s c 2 is the vertically integrated 
gas pressure, $ is the gravitational potential including 
the disk self-gravitational potential (if considered), n s is 
the viscous stress tensor for gas, respectively. The terms 
Ej n and F; n indicate the changes in the equations due to 
the infall model, where Ej n is the mass infall rate and 
Ei n is the shear force. We note that in our standard 
model we look at the effect of density enhancement only 
by matching velocities of infalling material and disk ma¬ 
terial (.F; n = 0). The shear term F- m in the momentum 
equation is added later in the model where the shear in 
between infalling material and disk material is considered 
(see below and Table []}. 

In the original work of [Ulrichl (|1976l ) and 
ICassen fe Moosmanl l|1981f l. the collapse of an isother- 

4 http://www.eso.org/public/news/esol436 


mal, spherically symmetric, and uniformly rotating 
cloud was studied. The infalling material follows 
parabolic orbits, arriving at the disk surface with 
different radial and azimuthal velocities from those 
of the disk material. Therefore, shear force exists 
which can be written as Fri n = 'F‘i n (vR,in — t’it.disk) and 
~ ^in(W,in ^.disk)- Here, = (GiU*/R) ^ 
and Vr/, : i n = (GM*/i? c ) -1 / 2 are the velocities of the 
infalling m ateri al with R c being the centrifugal radius 
(jCassen fe Moosmanl Il98lh . and u^uisk and v^^isk are 
the velocities of the disk, respectively. 

We consider tw o m odificat ions o f the infall mod el intro¬ 
duced in lUlrichl (|1976l) and ICassen fc Moosmanl (|1981l) . 
First, we simplify the model in a such way that the ra¬ 
dial and azimuthal velocities of the infalling material 
match those of disk material and thus no shear force 
exists (hereafter UCM m odel). On e notab le feature_of 
the infall model of lUlrichl (119761) and ICassen fc Moosmanl 
l|1981f) is that because of the solid-body rotation the in¬ 
falling material near the rotational axis has less angular 
momentum and thus falls at small radius while the in¬ 
falling material far from the rotational axis have more 
angular momentum and falls at large radius. This will 
concentrate infalling material near the outer edge of the 
infall (i.e. centrifugal radius). In order to alleviate the 
relatively strong density enhancement around the cen¬ 
trifugal radius, we consider another modification (here¬ 
after MUCM model). In the MUCM model, the radial 
infall pattern is modified in a such way that mass flux 
per unit distance is constant over radius in order to avoid 
the singularity in density of the infalling material at the 
centrif ugal radius (see below and Figure 1 of lBae et al.l 
l2013al for comparison between UCM and MUCM model). 

The mass infall rate of the UCM model is 

o-ir /2 “< 3 > 

and 

E in (R) = 0 if R > R c , (4) 

where M- ln = 0.975c;) c /G is the constant total infall mass 
rate at a given cloud isothermal sound speed c sc for the 
singular sphere solution (lShulll977l) . 

In the MUCM model, the mass infall rate is smoothed 
as 

E in (R) = A ^ r ‘ if R < R c (5) 

v ' 2nR c R ~ w 

and 

E in (R) = 0 if R > R c . (6) 

2.2. Dust Component 

In order to investigate the dust response to RWI- 
generate d ga s structures, we use the FARGO code 
(|Masse t 2000) that is modified to deal with two fluids as 
introduced in iZhu et all (I2012D . We treat the dust com¬ 
ponent as an inviscid, pressureless fluid and dust simply 
feels the drag force in addition to the central stellar po¬ 
tential. Note that the dust component in our calculations 
evolves passively and does not affect gas evolution and 
therefore the RWI. Dust feedback may become im por- 
tant at least locally in vortices (e.g. iFu et al.l I2014T) but 
the effect is not considered in our study. 
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TABLE 1 

Model Parameters 


Models 

Case Name 

a 

Rc 

(AU) 

Shear Terms 

Numerical Resolution 
(Nr X N$) 

Infall Model 

Self-gravity 

Standard Model (N3.lt 

s 

10- 4 

25, fixed 

N 

512 x 1024 

UCM 

N 

Numerical Resolution (N3.2t 

NR256 

10 -4 

25, fixed 

N 

256 x 512 

UCM 

N 


NR1024 

10" 4 

25, fixed 

N 

1024 x 2048 

UCM 

N 


NR2048 

10 -4 

25, fixed 

N 

2048 x 4096 

UCM 

N 

Viscosity (N3.31 

V2 

10-^ 

25, fixed 

N 

512 x 1024 

UCM 

N 


V3 

10 -3 

25, fixed 

N 

512 x 1024 

UCM 

N 


V5 

10" 5 

25, fixed 

N 

512 X 1024 

UCM 

N 

Increasing R c (N3.4t 

IRC2 

ur 4 

25, linearly increasing 

2 AU per 1000 years 

N 

512 X 1024 

UCM 

N 


IRC5 

icr 4 

25, linearly increasing 

5 AU per 1000 years 

N 

512 X 1024 

UCM 

N 

Shear Terms (93.511 

SH 

icr 4 

25, fixed 

Y 

512 X 1024 

UCM 

N 

Infall Model (93.611 

MUCM 

ur 4 

25, fixed 

N 

512 X 1024 

MUCM 

N 

Self-gravity (93.711 

SG 

ur 4 

25, fixed 

N 

512 X 1024 

UCM 

Y 


The drag terms are added in an additional source step 


as 


and 


dvR,d _ VR,d — VR,g 
dt t s 


dv^ d _ 

dt t s 


where Vd and v g denote dust and gas velocities, and t s is 
the dust stopping time. With the initial setup explained 
below in 1 )2.31 the mean free path of gas molecules is 
A = 3.8 (A/1 AU) 4 / 9 exp(A/A c ) cm so in our simula¬ 
tion domain the dust p articles smaller than ~ 1 m are 
in th e Epstein regime (Whippli Il972t iWeidenschillind 
2973). Therefore, the dust stopping time can be written 
as 


t s 


Pp s 


(9) 


PgVT 


where p p is the dust particle density, s is the dust par¬ 
ticle radius, p g is the gas density, and vt = -\/8/7rc s is 
the m ean thermal velocity w ith c s being the gas sound 
speed (jTakeuchi fc Linl [20021 1 . The dust particle density 
is assumed to be p p = 1 g cm -3 in this study. With 
the two-dimensional approach the gas mass density is 
p g = 'E g /y/2nH and thus the dust stopping time can 
also be written as 


, _ TtPpS 

s ~ 2 


( 10 ) 


where Q is the Keplerian angular velocity. The stokes 
number, or the nondimensional stopping time, T s is then 


T s = t s £l. 


( 11 ) 


In this work, we consider two different dust sizes: small 
dust particles that have T s ~ 0.1 and large dust particles 
that have T s ~ 1, under the initial conditions explained 
in the next section. In physical size, the former and latter 
corresponds to 1 cm and 10 cm, respectively. 


2.3. Initial Conditions 

We begin with a 0.5 M 0 central protostar and a sur¬ 
rounding disk having an initial gas density distribution 
of E(A) = 1000 (A/AU) _1 exp(— R/R c ) g cm -2 . Here, 


R c is the centrifugal radius which is fixed to 25 AU dur¬ 
ing the calculations unless otherwise stated. We assume 
a constant infall rate of 3.0 x 1CU 6 M^yr -1 which c or- 
responds to a singular isothermal collapse ( Shulll977 ) of 
a 16 K protostellar cloud. The disk mass is 0.014 M 0 
initially and is ~ 0.11 — 0.13 M 0 at the end of calcu¬ 
lations, depending on model parameters. We adopt a 
fixed radial temperature distribution (T rx A -1 / 2 ) that 
corresponds to a ratio of disk scale height to radius 
H/R = 0.05 (A/1 AU) 0 25 . Initial random perturbation 
is applied to surface density at the level of 10~ 4 . 

We use inner and outer boundaries of 5 AU and 100 AU 
and adopt 512 logarithmically spaced radial grid-cells 
and 1024 linearly spaced azimuthal grid-cells. With this 
choice, AA/A is constant to 0.006 and grid-cells have 
comparable radial and azimuthal size at all radii. In 
i )3.21 we discuss the effect of numerical resolution and 
show the results converge resolutions at 512 x 1024 and 
beyond. 

We use a = 10” 4 for our fiducial viscosity parame¬ 
ter. Using a relatively low disk viscosity is motivated by 
the fact that stellar X-rays and FUV photons are likely 
to be blocked by infalling material during the protostel¬ 
lar phase so that ionization level at disk surface layers 
remains low, limiting mass transport through the mag- 
netorotational instability. The effect of the viscosity pa¬ 
rameter is tested in i )3.3l Model parameters are summa¬ 
rized in Table 21 

2.4. Boundary Conditions 

In order to prevent unphysical, rapid depletion of ma¬ 
terial at the inner boun dary, we implemen t a ve locity 
limiter as intro duced in IPierens k, Nelsonl (120081) and 
iZhu et all (|2012fb To briefly summarize, we limit the ra¬ 
dial velocities of the gas component at the inner bound¬ 
ary to be no more than some factors of the viscous radial 
velocity in a steady state, VRy ls : vr >3 < /3vr :V i S . Here, 
VR,v is = 3iZj n /27?|j^ where , and R[n are the viscosity 

and radius at the inner boundary. Based on a set of ex¬ 
periments, we find that /3 = 3 is most suitable for this 
study. 

A similar appr oach is implemented for the dust com¬ 
ponent, following IZhu et abl (|2012l ). In case of dust, we 
limit radial velocity to be no more than three times of the 
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Fig. 1.— Radial distributions of azimuthally-averaged (a) gas 
surface density £ ff , (b) gas azimuthal velocity v^ } g, (c) epicyclic 
frequency and (d) vortensity r}/rj(R c ) at the launching 

of the RWI (t = 14 orbital times at R c ; see Figure [2} for the 
standard model. Dashed curves show the initial distributions. The 
centrifugal radius R c is indicated with vertical dotted lines. In 
panel (b), two dotted-dashed lines show slopes of —1 and —0.5, 
respectively. We note that the density bump generates a steep 
azimuthal velocity gradient around R c in order to maintain the 
radial force balance, which results in the k, 2 and 77 minima. 


Fig. 2.— The maximum value of non-axisymmetric gas density 
perturbation SH g /(Y! l g) is shown as a function of time. The insta¬ 
bility triggers at 14 T or b and grows exponentially during the ‘linear’ 
phase (14 T or b < t < 30 T or b). Then the instability saturates and 
turns into the nonlinear regime thereafter. 

P/E 7 is the entropy. Given that we assume a locally 
isothermal temperature profile which is unchanged over 
time, the above equation can be written as 


dust drift speed in a viscous disk 'C/?, drift; where Adrift 
is defined as 


(3iy in /2R in )T s 1 +r/v K 

Adrift- Ts+T _ 1 


( 12 ) 


Here, 77 is the ratio between the pressure gradient and 
gravitational force defined as 77 = —(RCl 2 p g )~ 1 dP/dR, 
and vk is the Keplerian speed. 

At the outer boundary, standard open boundary condi¬ 
tions are implemented for both gas and dust components. 


3. RESULTS 

3.1. Standard Model: Results With the UCM Model 

We start with our standard model, showing details of 
the launching, growth, and saturation of the RWI and 
dust’s response. Then, in the following subsections we 
present how varying disk conditions and infall pattern 
affect the RWI activity. 


K 2 1 

2 E n cf' 


(14) 


For convenience we call the quantity 77 as vortensity here¬ 
after. 

Figure |Tj shows radial distributions of azimuthally- 
averaged gas surface density, azimuthal velocity, epicyclic 
frequency, and vortensity at the time of the launching of 
the RWI. The density gradient and thus the pressure gra¬ 
dient becomes steeper on the outer part of the density 
bump as infall adds mass onto the disk. The steep pres¬ 
sure gradient then generates a steeper azimuthal velocity 
profile in order to maintain the radial force balance. The 
resulting azimuthal velocity slope deviates from —0.5, 
which is for the Keplerian rotation, and approaches to 
— 1. We note that the —1 slope is important in trigger¬ 
ing the RWI since n 2 and 77 change sign at oc i? _1 . 
This can be simply shown by substituting n 2 into Equa¬ 
tion (fill) : 


3.1.1. Launching of the RWI 

In two-dimensions, one can analytically show that the 
RWI can be triggered at a radial minimum of the gener¬ 
alized vortensity 77 , where 

« 2 1 /.of 

V 2EU S 2 h ^ ) 

(|Lovelace et all 119991 : ILi et al.1 1200011 . Here, n — 
[R~ 3 d(R 4 V, 2 )/dR] 1 ^ 2 is the epicyclic frequency, E is the 
surface density, Q is the rotational frequency, and S = 


1 = 


111 d 
2EUcf R^dR 


(R 2 vD- 


(15) 


In Keplerian disks where v^ oc i? -1 / 2 , 77 is always posi¬ 
tive. However, 77 becomes zero if v</, oc R~ l and is nega¬ 
tive if has a steeper slope than — 1 . 

Figure [2] presents the maximum value of non- 
axisymmetric density perturbation <5E g /(E g ) as a func¬ 
tion of time, where <5E g = E g — (E g ) and the brackets de¬ 
note the azimuthal average. Initially before the RWI trig¬ 
gers, the maximum perturbed density maintains ~ 10~ 4 


























ROSSBY WAVE INSTABILITY BY PROTOSTELLAR INFALL 


5 


< 


>- 







- 2.0 


- 1.0 


EH xio -4 H 

2.0 - 0.65 




- 0.40 


- 0.15 


0.10 


EH 

0.35 




- 0.65 - 0.40 - 0.15 


0.10 


EH 

0.35 


-50 


-50 


-50 


-100 
-100 


-50 


0 

X [AU] 


100 


-100 
-100 


-50 


100 


-100 
-100 


-50 


100 


Fig. 3.— Perturbed gas density distributions 5E g /(E g ) of the standard model at the launching of the instability, at the saturation, and 
at the end of the simulation (from left to right). These correspond to 14, 30, and 226 local orbital times at R c or t = 2.5, 5.3, and 40 X 10 3 
years. Note that the scale for the leftmost panel is different from the other two. 


which corresponds to the initial random component. The 
vortensity minimum develops gradually until the insta¬ 
bility triggers at t ~ 14 T or b, where T or b hereafter denotes 
the local orbital time at the centrifugal radius. 

3.1.2. Growth and Saturation of the RWI 

After the instability is triggered, it grows in the ‘lin¬ 
ear’ regime where the growth of the perturbed density 
is well described by the linearized continuity equation 
d(8Yi)/dt = —iCSSYi, where Co = lo + i'y is a complex 
frequency, w is the real mode frequency, and 7 is the 
growth rate. This is clearly seen in Figure [2] the 
instability grows exponentially during the linear phase 
(14 T orb < t < 30 T orb ). The growth timescale T gr 0 wth, 
during which time the maximum perturbed density in¬ 
creases by a factor of e, is 1.8 T orb . After spending 16 T orb 
in the linear phase, the instability saturates and enters 
the nonlinear regime. 

Figure [3] shows the spatial distributions of the per¬ 
turbed gas density at the launching of the instability 
{t = 14 T orb ), at its saturation (t = 30 T orb ), and at 
the end of the simulation (t = 226 T orb ). Initially, the 
m = 4 mode is dominant but it quickly merges to to = 3 
mode during the linear phase, and eventually merges to 
to = 1 mode which is maintained until the end of the 
calculation. 

3.1.3. Vortex Formation and Dust Response 

The RWI accompanies vortex formation as the insta¬ 
bility enters the nonlinear regime. Figure [4] illustrates 
gas vorticity, gas surface density, and surface densities of 
1 cm and 10 cm dust particles on the cj> — R coordinates. 
The (negative) vorticity minimum grows at R c as infall 
proceeds. At the time the RWI initiates (t = 14 T orb ), 
the non-axisymmetric features are still too small to be 
seen. As the instability grows, the radial vorticity mini¬ 
mum starts to develop structure and finally breaks into 
vortices when the instability enters the nonlinear regime 
(t = 30 Tc-t,). Vortex formation during nonlinear evo¬ 
lution of the RWI is in good agreem ent with previous 
hydrodynamic simulations fe.g. iLi et al.ll200lD . The im¬ 
portant feature is that the vortices have a local vorticity 
minimum at its center - the vortices are anticyclones. 


TABLE 2 
Results 


Case Name 

RWI 

K 2 /n|. a 

Yiaunch b 

T sa t b 

rp b 

- 1 growth 

Vortex 



(To rb) 

(Torb ) 

(T 0 rb) 

Formation 

s 

Y 

0.16 

14 

30 

1.8 

Y 

NR256 

Y 

0.24 

13 

29 

2.1 

Y 

NR1024 

Y 

0.15 

14 

30 

1.7 

Y 

NR2048 

Y 

0.15 

16 

32 

1.7 

Y 

V2 

Y c 

0.62 

22 

- 

2.2 

N 

V3 

Y 

0.25 

18 

36 

2.1 

Y 

V5 

Y 

0.15 

14 

30 

1.8 

Y 

IRC2 

Y 

0.28 

22 

33 

2.2 

Y 

IRC5 

Y c 

0.60 

24 

- 

5.3 

N 

SH 

Y 

0.03 

8 

19 

1.2 

Y 

MU CM 

Y 

0.35 

54 

110 

5.5 

Y 

SG 

Y 

0.15 

14 

32 

1.9 

Y d 


a The k 2 /£l 2 K values are measured at the launching of the RWI. 
k The times correspond to the launching of the RWI (Ti aU nch), the saturation of 
the RWI (T S at)i and the exponential growing timescale (T growt h)- 
c The instability is triggered but stays only in the linear regime and withers away 
before it enters the nonlinear phase. 

^ Multiple vortices form as the RWI enters the nonlinear regime, but they later 
dissipate as the disk becomes gravitationally unstable. 

The vortices wander in azimuth, merge together, and 
form a single vortex within ~ 20 T orb from their for¬ 
mation. The merged vortex survives until the end of the 
calculation and we did not follow its evolution thereafter. 

Before the vortices form, dust particles are concen¬ 
trated at the vorticity minimum while the small and the 
large dust behave somewhat differently. Initially when 
the instability grows in the linear phase, small dust con¬ 
centrate around the gas pressure maximum. On the other 
hand, the large particles have the greatest inward drift 
velocity at the centrifugal radius (T s ~ 1 at i? c ), so the 
density drops at the radius. The large dust distribution 
shows a steep density change around R c (see Figure [4] 
at t = 14 and 26 T orb ). This is because infall adds gas 
inside R c and thus the dust stopping time correspond¬ 
ingly decreases significantly at the region, slowing down 
the large dust migrating inward. 

After the vortices form, they efficiently trap dust parti¬ 
cles so that the dust distributions show stronger asymme¬ 
try than the gas distribution. While the small particles 
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Fig. 4.— Snapshots of gas vorticity V x v g , gas density, 1 cm particle density, 10 cm particle density at selected times. Times are 
presented on the left side in units of T or b. The Keplerian component in the azimuthal velocity is subtracted when calculating the vorticity. 
Densities are in cgs units and displayed in the logarithmic scale. 


are well coupled to gas and follow the gas vorticity dis¬ 
tribution well, the large particles show highly perturbed 
behavior. The accumulated cores of large dust are often 
offset from the vortex cores. This is more significant in 
earlier times when the vortices merge (see Figure |4] at 
t = 38 and 40 T or b for example). Also, the large dust 
particles show more dramatic concentration than small 
dust particles. We will discuss dust trapping and its im¬ 
plication later in 1 14.11 


3.2. Effect of Numerical Resolution 

As mentioned in EH the UCM model has a singu¬ 
larity (infinite mass infall rate per unit area) at the 
centrifugal radius while the total infall rate is finite. 
With finite grids this can cause non-convergent behav¬ 
ior at different numerical resolutions. Thus, we need to 
show that our results are not dependent on the numeri¬ 
cal resolution before we go to further analysis. We test 
with four different numerical resolutions of {Nr, Nff) = 
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Fig. 5.— (a) Radial distribution of the azimuthally-averaged r]/r](Rc) at t = 14 T or b- (b) Maximum perturbed gas density <5£ g /(E g ) 
as a function of time, (c) The Reynolds stress as a function of time. In panel (b) and (c), the plots are shifted vertically by 1 and 0.002, 
respectively, for better view. We note that the oscillating feature of QR ey is due to the epicyclic motion of gas inside vortices against the 
geometric center of the vortices. Since it requires a lot of computational resources the highest resolution run is conducted only for 80 T or b. 
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Fig. 6.— Radial distribution of azimuthally-averaged r]/r](Rc) at 
t = 14 T or b with a = 10“ 2 , 10“ 3 , 10 —4 (standard model), and 
10 -5 . The vortensity minimum is very shallow and broad with 
a = 10~ 2 . In this case, the RWI grows only in the linear regime 
and does not further develop to the nonlinear regime. 


Fig. 7.— Radial distributions of the azimuthally-averaged 
rj/ri(Rc) with fixed R c at 25 AU (standard model) and linearly 
increasing R c at rates of 2 AU and 5 AU per 1000 years. Although 
the vortensity minimum is seen in all models, the RWI grows only 
in the linear regime when the centrifugal radius increases 5 AU per 
1000 years. 


(256, 512), (512,1024), (1024, 2048), and (2048,4096). 

We use three different diagnostics to check the numeri¬ 
cal convergence. First, we check the position, shape, and 
depth of the vortensity minimum since it is crucial for 
triggering the RWI. Second, we use the maximum value 
of the perturbed gas density. For the last, we calculate 
the Reynolds stress which is defined as 


_ / T,g5VR t g5v < l >t g dS 

aR ^~ JS^.dS ’ 


(16) 


where the integration is done in a volume-averaging man¬ 
ner, and Sv Rt g = v R g - (v R , g ) and <5- ( v 0 , s ). 
The first two quantities check the convergence of local 
features whereas cxR e y checks the numerical convergence 


in a globally-averaged sense. 

The three diagnostics are presented in Figure [5] All 
diagnostics converge well at 512 x 1024 and beyond. With 
256 x 512 grid cells, the region around the R c is not very 
well resolved so that the infall rate seems underestimated 
at the region. The vortensity minimum with 256 x 512 
grids is significantly shallower and broader than others 
(see also Table [2] for the minimum h 2 /VL 2 k values), and 
thus the instability growth and the Reynolds stress are 
less significant. 


3.3. Effect of Viscosity 

Since larger gas viscosity more efficiently spreads out 
the density enhancement, it is possible to broaden the 
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Fig. 8. — Same as Figure [l] but at t = 8 T or b for the SH model 
(shear terms are included). We emphasize that the vortensity min¬ 
imum is much sharper and narrower compared to the standard 
model (see Figure [TJl). 


vortensity minimum and limit the RWI growth. We 
tested different viscosity parameters to see the effect of 
disk viscosity; a = 10 -2 ,10~ 3 ,10 -4 , and 10 -5 . As ex¬ 
pected, Figure [S] shows that the vortensity minimum be¬ 
comes broader as a increases. Especially with a = 10” 2 , 
the minimum is very shallow and broad. In this case, the 
RWI triggers but it only stays in the linear regime and 
does not turn into the nonlinear regime - we therefore 
do not observe any vortices forming and concentration of 
dust particles. 

It has been empirically shown that the width of the 
vortensity minim um has to be < 2 H in order for th e 
RWI to develop (ILvra et all 120091 IReealv et all 120121) . 
If we adopt this criterion, the RWI growth can be lim¬ 
ited if the local viscous timescale at R C: in which time 
the disk can viscously spread out the bump by ~ H, 
is shorter than the RWI growing timescale. The lo¬ 
cal viscous timescale estimated at R c is t v = H 2 /v ~ 
(0.11 R) 2 /v ~ 90/a years. With a = 10 —2 , the local vis¬ 
cous timescale is ~ 9000 years (or ~ 51 T or b) which is 
approximately the timescale of the linearly growing part 
of the RWI. Thus, it makes sense that the instability 
with a = 10 -2 does not further grow and turn into the 
nonlinear regime. On the other hand, viscosity does not 
affect triggering of the RWI when the viscous timescale 
is much longer than the RWI growing timescale. 

3.4. Effect of Linearly Increasing Centrifugal Radius 
With Time 

In protoplanetary disks, the centrifugal radius will gen¬ 
erally not be constant with time as material with even 
greater angular momentum falls in; this will spread the 
vortensity minimum. Instead of implementing a more 
self-cons istent time ev olution of the centrifugal radius 
(e.g. iCassen fc Moosmanlll98il h we mimic the effect by 


simply linearly increasing R c as a function of time. We 
test with two different rates: R c increases at 2 AU and 
5 AU per 1000 years. 

In order to see how the R c increasing rates translate 
to initial core rotation, let us assume that the proto- 
stellar system we consider in this study evolves from 
the initial Bonn er-Ebert sphere- like two-component den¬ 
sity profile (see lZhu et alll20Tcil) . and the inner flat core 
has collapsed to a 0.1 M 0 central protostar and the 
rest of the cloud collapses as a rotating singular isother¬ 
mal sphere. If we further assume that the disk mass 
inside of our inner boundary (R ln = 5 AU) is negligi¬ 
ble, our model is ~ 1.3 x 10 5 years past from the ini¬ 
tial core collapse. Then, the increasing rates of R c we 
adopt here correspond to the instantaneous R c increas¬ 
ing rate of uniformly rotating cores at fl c = 5.3 x 10 -14 
and 8.3 x 10“ 14 rad s _1 . The rotation frequencies are 
20 and 30 % of the breakup angular frequency at the 
outer edge of an 1 M 0 cloud. It is worth to point out 
that the rotation frequencies are large com pared to the 
median value of 3 % inferred bv IBae et all (I2013bfl : the 
value was required to reproduce the observed circumstel- 
lar disk frequencies as a function of age, assuming disk 
dispersal by photoevaporation. 

Figure [7] shows radial distributions of vortensity mini¬ 
mum for the models. The RWI excites and generates vor¬ 
tices with small increasing rate of 2 AU per 1000 years. 
However, if R c increases fast at 5 AU per 1000 years the 
RWI triggers but stays only in the linear regime. The 
instability with the large rate eventually withers away 
and no vortices form in this case. Although a more self- 
consistently evolved model is require to conclude, in very 
rapidly-rotating systems the centrifugal radius may move 
outward so fast that the RWI may not have chance to 
form vortices. 


3.5. Effect of Shear Terms 

In actual protostellar systems, the infalling material 
must have different specific angular momentum or az¬ 
imuthal velocity than the disk material when it lands on 
the disk. Otherwise, the mass would not fall in. We 
test the effect of shear by adopting velocity fields of in- 
falli ng m aterial that foll ows parabo li c orb its. Follow¬ 
ing lUlrichl (|1976T ) and ICassen fc Moosmanl (|1981l h the 
infalling material has radial and azimuthal velocities of 
VR, in = -(GM./fl)- 1 / 2 and v^ n = (GM*/i? c ) -1 / 2 . So 
at R = R c both radial and azimuthal velocities of in¬ 
falling material are the same as the Keplerian velocity at 
the radius. We note that the radial velocity of infalling 
material is extremely large when compared to the accre¬ 
tion velocity in the steady-state a disk. In an a disk, the 
accretion rate can be described as M acc cs SirnT, since 
we are interested in the region far from the stellar sur¬ 
face. Then, if we relate it to M acc = — 27rI?£?;R i disk we 
get VR^isk = -(3/2 )ac s (H/R). Thus, i>R,di s k/FR,in = 
(3/2 )a{H/R) 2 < 1. 

Due to the large inward radial velocity of the infalling 
material, the disk develops a very sharp density jump 
at the centrifugal radius. Therefore, compared to the 
standard model where no shear terms are included, the 
vortensity minimum is much sharper and narrower as 
seen in Figure [8] In this model the shear from infall 
rapidly builds the vortensity minimum, whereas in the 
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Fig. 9.— Perturbed gas density <5£ g /(S g ) distributions at the launching of the instability, at the saturation, and at the end of the 
simulation (from left to right) for the SH model (shear terms are included). These correspond to t = 8,19, and 226 T or n or 1.4,3.4, and 
40 X 10 3 years. Note that the scale for the leftmost panel differs from the other two. 


standard model the disk has to wait until enough mass is 
added to build the density bump. Thus, no evident den¬ 
sity bump is developed around R c at the time of RWI ini¬ 
tiation. Another notable feature is that the RWI initially 
triggers with an extremely high order mode of to = 9 as 
seen in Figure |9l at which mode the linear growth rate is 
the highest in this model. The perturbed density peaks 
rotate at slightly different velocities and thus one catches 
another as time goes. They eventually merge to m = 1 
mode at i ~ 100 T 0 rb- 

3.6. Effect of Infall Profile: With the Modified UCM 

Model 

The UCM model, as pointed out earlier, adds a large 
fraction of infalling material near the centrifugal radius. 
In order to see if the RWI can be limited by more gentle 
infall pattern, we match the radial infall pattern to the 
initial disk surface density distribution (E g , E; n oc i?” 1 ). 
Even with the smoothed infall profile we find that the 
RWI excites. Figure [TO] presents radial distributions of 
azimuthally-averaged gas surface density, azimuthal ve¬ 
locity, epicyclic frequency, and vortensity at the time of 
the launching of the RWI (t = 54 T or b). The vortensity 
minimum is shallow and broad, but since infall keeps 
adding material at the same radius the RWI triggers. 
The density around R c increases steeply but does not 
develop bumpy structures. We note that the instabil¬ 
ity very slowly grows, having three times longer T growt h 
compared to the standard model (see Table [2J. 

3.7. Effect of Self-gravity 

In our standard run, the azimuthally-averaged Toomre 
Q parameter around R c is < 1 at T > 200 T or b. At 
this point (or even earlier), we expect disk self-gravity 
becomes important and may alter the later disk evo¬ 
lution including possible activation of gravitational in¬ 
stability (GI). We thus conducted a calculation with 
disk se lf- grav ity inclu ded , using FARGO-ADSG code 
(iBaruteau & Massetll2008f ). 

Before looking at numerical results, one may predict 
the role of self-gravity in triggering the RWI through a 
simple back-of-the-envelope calculation using the connec¬ 
tion between the vortensity 77 and the Toomre Q param¬ 





Fig. 10.— Same as Figure |T] but for the MUCM model. 


eter, where the vortensity is again 

K 2 t 

V ~ 2UE cf ’ 

and the Toomre Q parameter is 


Q = 


KC S 

7rGE' 


(17) 


(18) 


At a given radii, under the locally isothermal assumption, 
77 oc k 2 /E and Q oc k/ E. Remember that both RWI 
and GI acts in a way to redistribute mass in the disk 
so the disk stabilizes against the instabilities. In other 
words, the instabilities broaden the density enhancement 
and increase the epicyclic frequency close to Keplerian 
speeds. If a disk is under the circumstance that RWI 
and GI competes, this process will increase 77 faster than 
Q so that the RWI will stabilize first. Therefore, one can 
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Fig. 11.— Perturbed gas density <5S g /(E g ) distributions at /. = 30,44, and 160 T or n (from left to right) for the standard run (upper 
panels) and the self-gravity run (lower panels). The azimuthally-averaged Toomre Q parameters of the self-gravity run at the three epochs 
are 5.3, 3.2, and 1.3. 


expect RWI will still trigger while a disk is gravitationally 
stable (Q> 1) but as the disk becomes gravitationally 
unstable (Q ~ 1) GI can eventually take part in, in which 
case the RWI can be quenched. 

Figure [FT] depicts perturbed density distributions at 
three selected times (t = 30,44, and 160 T or b), together 
with the distributions of the standard run for compar¬ 
ison. During the period the RWI linearly evolves, the 
azimuthally-averaged Toomre Q parameter in the disk is 
> 5; the disk is gravitationally stable. The RWI thus 
triggers as we predict above and vortices form. Vortices 
form as the RWI enters the non-linear regime and merge 
together over time. However, in the presence of self¬ 
gravity the merging tends to be impeded and results in 
m = 2 mode, instead of a single vortex. At t = 44 T or b, 
azimuthally-averaged the Toomre Q parameter is 3.2 at 
the radial density bump while Q is locally as small as 
2.4 at the core of the vortices. As the disk further ob¬ 
tains material and becomes gravitationally unstable, the 
vortices dissipate and the GI eventually triggers gener¬ 
ating strong m = 2 trailing spiral arms. The evolu¬ 
tion of vortices under the influence of self-gravity seen 
in our calculation agrees well with previous two- and 
three-dimensional simulations showing that vortex merg¬ 
ing can be delayed in weakly self-gravitating disks and 
global spiral waves develop inste ad of vortex formation if 
self-gravi ty is sufficiently strong (|Lin fe Papaloizoull20lTI : 

E^EqH). 
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Fig. 12.— Time evolution of the minimum gas to dust ratio for 
the standard model. The ratio drops significantly as the instability 
enters the nonlinear regime and vortices form. Because the large 
dust at the outer disk are depleted due to rapid radial drift and 
there are no more supply from the region, the ratio for the 10 cm 
particles after ~ 120 T or b increases. 


4. DISCUSSION 
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Fig. 13.— Surface density distributions of (left) gas, (middle) 1 cm dust particles, and (right) 10 cm dust particles at the end of the 
simulation for the standard model. Densities are in cgs units and displayed in the logarithmic scale. 


4.1. Dust Trapping in Vortices and Observational 
Implications 

Our two-fluid approach shows that vortices formed by 
the RWI efficiently trap dust particles of appropriate size. 
Figure |T2] shows the time evolution of the minimum gas 
to dust ratio for the standard model. Although the de¬ 
tailed evolution of the ratio for the small (1 cm) and 
large (10 cm) dust particles are different, in overall, the 
ratios tend to decrease by more than a factor of few. The 
small dust particles are well coupled to gas as we see in 
33.1.31 and therefore the decrease of the ratio is relatively 
smooth and only moderate. The ratio has a minimum of 
~ 20 : 1. For the large dust particles, the ratio slightly 
increases initially because they migrate inward rapidly 
(remember T s ~ 1 at R c ), while gas surface density in¬ 
terior of R c increases due to the infall. As vortices form, 
however, they efficiently trap the large dust particles and 
the gas to dust ratio drops to ~ 2.7 : 1 at t ~ 120 T or b- 
The ratio then increases after ~ 120 T or b but this is be¬ 
cause the large dust at the outer disk are depleted due 
to rapid radial drift and there are no more supply from 
the region. 

Figure H31 displays surface density distributions of gas, 
small dust, and large dust at the end of the standard 
run. Clear dust concentration in the vortex is shown 
for both small and large dust particles. Another no¬ 
table feature is that the inner disk is depleted of dust 
because of the dust trapping by the vortex. Again, this 
effect is more prominent for the larger dust particles. 
As seen in the figure, existence of vortices can not only 
form azimuthal asymmetry of dust that can be directly 
observed via interferometric observations but also build 
radial structures that can be inferred by single-dish ob¬ 
servations at multiple wavelengths. The deficit of dust 
at inner disk regions can affect the spectral energy distri¬ 
butions of the disk and su ch d isks can be interpreted as 
transitional disks dCalvet et alJl2005h . The disks still can 
have enough gas at the inner disk (see Figure [13]) that 
can probably maintain high accretion rate observed in 
some transitional di s ks (e.g. lEspaillat et al.l 1201(1120111 : 
lAndrews et al.ll201lt iKim et al.l 120131 ). 

The significant dust concentration in vortices suggests 
that the vortices can provide favorable conditions for the 
planet and/or planetesimal formation. In terms of vor¬ 


tex formation via the RWI, it would prefer higher infall 
rate and lower disk mass although a thorough parameter 
study is desired in the future regarding infall rate and 
disk mass. This means that, with the example of HL 
Tau suggesting that planet formation can start early even 
during the infall phase, RWI might be one mechanism 
to accelerate planet formation. On the other hand, the 
longevity of vortices has to be fu rther tested since vor - 
tices may not survive forever (e.g. lMeheut et a l. 2012b), 
although dust c an still remain concen trated after gas vor¬ 
tices disappear (jBirnstiel et al.ll2f)13l ). 

4.2. Angular Momentum Transport 

iLi et all ( 200 1) has shown that RWI-driven anticy¬ 
clones and trailing spiral waves are responsible for out¬ 
ward angular momentum transport. In Figure [14] we 
display spatial distribution of the perturbed gas den¬ 
sity along with the radial distribution of azimuthally- 
averaged Reynolds stress. At the time the instability 
saturates, the dominant m = 3 mode is clearly seen at 
all radii: the trailing spiral waves propagate both inte¬ 
rior and exterior of R c . The measured Reynolds stress 
has the maximum value of ~ 0.015 around R c . Inside 
of R c the Reynolds stress is measured to ~ 10 -3 while 
at R > 60 AU the stress is in the range of aR ey ~ 
10 -4 —10~ 6 . At later time, the vortices merge and the hy¬ 
drodynamic turbulence around the vortex becomes less 
significant; QR e y ~ 10 -3 . Instead, the trailing one-armed 
spiral wave generate stronger turbulence very broadly in 
the disk that corresponds to dR ey ~ 10~ 4 — 3 x 10 -3 . 
The measured Reynolds stress driven by the RWI and 
the subsequent vortex formation is as strong as the ones 
induced by other instabilities in protoplanetary disks, 
such as gravitational instability (e.g, iBae et al.ll 2014) and 
magnetorotational instability ('e.g. lStone et al.l!l996D . 

As a result of exchange of angular momentum between 
the vortex and the surroundin g disk mater ial, t he vor - 
tex is subject to migration (jPaardekooper et all (2010 ). 
More specifically, the vortex loses angular momentum 
and migrates inward if the trailing waves emitted by the 
vortex are stron ger at the outer disk reg ion than the in¬ 
ner disk region (jPaardekooper et al.ll20ldh . We measure 
the migration rate starting from the time vortices merge 
together in a single, stable vortex, which is t ~ 80 T or b 
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Fig. 14. — (Top) Perturbed gas density <5E g /(E g ) distributions on the R — (f> coordinates. (Bottom) Azimuthally-averaged Reynolds stress 
as a function of radius. The snapshots are taken from the standard model at t = 30 T ot u (left) and 226 T or b (right). 


in the standard model for example (see Figure [fj. How¬ 
ever, we do not see the vortex migrating until the end 
of the simulation in our models. It has to be noted 
though that this can be due to finite numerical resolu¬ 
tion. Since the radial grid size at R c is ~ 0.15 AU, in 
the standard run we are able to observe migration only 
if the rate is > 6 x 10~ 6 AU yr” 1 . We also do not ob¬ 
serve vortex migrating in the higher resolution run with 
1024 x 2048 grid cells, which reduces the migration rate 
to<3xlO _6 AU yr _1 . The fact that the vortex migrates 
only very slowly, if it does, is important because the vor¬ 
tex might prevent centimeter to meter-sized objects from 
rapidly migrating inward in a few x 10 3 — 10 4 years via 
the aerodynamic drag, retaining material for planet for¬ 
mation. 

Generally, it is known that v ortices mi gra t e to¬ 
ward high pressure regions (|Pa ardckoo i)er et all 120101 : 
iMeheut et al.l l2012bl) . If a disk around a vortex has 
constant pressure, competing effects that determine vor¬ 


tex migration cancel each other an d th u s the vortex is 
likely to be locked in its location dPaa rdekoo per et al.l 
120101 : iLvra fe Mac Lowll2012t IMeheut et al.ll2012bl) . In 
our models, we find that the radial pressure gradient 
around merged vortices is shallow, in many cases close to 
zero, and this is presumably why we do not observe vor¬ 
tices migrating. However, we caution that migration of 
vortices could depend on many other complications (e.g. 
iRichard et al.ll2013tlFaure et al.ll2015l) . which are beyond 
the scope of two-dimensional adiabatic calculations. 

4.3. Caveats and Future Work 

It is likely that actual patterns of protostellar infall will 
be much more complex, possibly having fila mentar y infall 
pattern, as inferred from observa tions (e.g. iTobin et all 
12010112011L l2012t I Yen et ap 1 20 141) and suggest ed by nu¬ 
merical simulations (e.g. iSeifried et al.1 1201511 . These 
more complex patterns of mass and angular momentum 
addition might yield differing structures and other in- 
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stabilities including RWI. Three-dimensional simulations 
are needed to take the next steps toward understanding 
the development of structure in protoplanetary disks. 

It is possible that vortices generated during the proto- 
stellar infall phase dissipate during the subsequent disk 
evolution, especially under the influence of disk self¬ 
gravity as we show in 1 )3.71 In our example, however, the 
disk becomes gravitationally unstable quickly because we 
add infalling material at the same disk regions over time, 
whereas infall is likely to occur at larger radii as time 
passes, due to the addition of m ater ial with higher angu¬ 
lar momentum le.g. iCassen fc Moosmanlll981 L Further 
studies of the longevity of vortices formed during the in¬ 
fall phase which incorporate a self-consistent and realistic 
evolutionary model are thus required. 

Some other limitations of the current work include lack 
of the three-dimensionality of the RWI and vortex struc¬ 
ture, proper thermodynamics which is especially impor¬ 
tant when self-gravity is considered, and more accurate 
dust physics such as dust growth and dust feedback. 


5. CONCLUSION 

We propose protostellar infall as a possible mechanism 
to trigger the RWI. Our work demonstrates that the RWI 
enables early vortex formation in protoplanetary disks, 
during the infall phase. This, along with the emerging 
observation evidences may suggest that planet/planetary 
core formation can start earlier than previously expected. 

By implementing infall models, we carry out two-fluid, 
two-dimensional global hydrodynamic simulations. Our 
results show that the RWI triggers at a radial minimum 
of the vortensity which is in agreement with previous 
works. In our model the vortensity minimum develops 
near the density enhancement at the outer edge of the 
mass landing on the disk (centrifugal radius). The key 
feature of triggering the RWI is the steep radial gradient 
of the azimuthal velocity close to R~ x instead of R~ 0 5 
for Keplerian rotation, which is induced by the local in¬ 
crease in density at the centrifugal radius. The instability 
initially grows in the linear regime where the growth is 
well described by the linear theory. The vortensity mini¬ 
mum keeps growing as infall proceeds, and the instability 
eventually saturates, followed by subsequent nonlinear 
evolution. 

We conduct a parameter study to investigate the RWI 
activity under a variety of disk conditions. The major 
findings are as follows: (1) the RWI triggers with disk 
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